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Abstract. We study the Bethe Ansatz/Ordinary Differential Equation (BA/ODE) 
correspondence for Bethe Ansatz equations that belong to a certain class of coupled, 
nonlinear, algebraic equations. Through this approach we numerically obtain the 
generalised Heine-Stieltjes and Van Vleck polynomials in the degenerate, two-level 
limit for four cases of exactly solvable Bardeen-Cooper-Schrieffer (BCS) pairing models. 
These are the s-wave pairing model, the p + zp-wave pairing model, the p + ip pairing 
model coupled to a bosonic molecular pair degree of freedom, and a newly introduced 
extended d + id-wave pairing model with additional interactions. The zeros of the 
generalised Heine-Stieltjes polynomials provide solutions of the corresponding Bethe 
Ansatz equations. We compare the roots of the ground states with curves obtained 
from the solution of a singular integral equation approximation, which allows for a 
characterisation of ground-state phases in these systems. Our techniques also permit 
for the computation of the roots of the excited states. These results illustrate how the 
BA/ODE correspondence can be used to provide new numerical methods to study a 
variety of integrable systems. 
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1. Introduction 

The correspondence between zeros of polynomials and models of physical systems has 
a very long history which goes back to early works of Stieltjes [TJ, Heine [2], Bocher [3], 
Van Vleck [4] and Polya [5] at the end of the 19th century in relation to electrostatic 
models. A concise summary can be found in [6]. The subject of zeros of polynomials [7], 
and specifically the sum rules of their zeros, was systematically studied for Hermite, 
Laguerre, Tchebycheff, Jacobi and Lame polynomials [HUH]. The study of the generalised 
Lame equation (a second order Fuchsian equation), its Heine-Stieltjes and Van Vleck 
polynomials, and their zeros, is still a very active area of research [TOl - TlT] . 

The relation between ordinary differential equations and integrable systems has 
been observed by several authors [T8H22] and in particular connection to a generalised 
Stieltjes problem has been discussed [231121] . Richardson's Bethe Ansatz solution [251126] 
for the s-wave pairing Bardeen-Cooper-Schrieffer (BCS) Hamiltonian has attracted 
considerable attention [2?14"4"T] . The correspondence between the Richardson's equations 
and the confluent Heun differential equation was recognised by Gaudin [27]. Numerical 
methods to solve Richardson's equations for finite number of particles were implemented 
later [32]. Here there may be convergence problems concerning critical points where the 
solution set of roots contains degeneracies [36] • The numerical methods of [321136] rely 
on a direct approach, or use an appropriate change of variables for the critical points, 
to solve the Bethe Ansatz equations. 

With regard to the difficulty of directly solving the Bethe Ansatz equations, very 
recently new numerical approaches based on the Bethe Ansatz/Ordinary Differential 
Equation (BA/ODE) correspondence were proposed [39, H2TTH] . They are based 
on linear second-order differential equations [4*31 14*1"] . or their corresponding Riccati 
equations "3"9"1I4*2] which are first-order nonlinear differential equations. The polynomials 
obtained by these methods correspond to extended Heine-Stieltjes polynomials [23J. 
Many advantages of these approaches were discussed in these works and it appears that 
such methods could be beneficial for studying integrable systems at their critical points. 
However, the applications of these new methods were limited to a certain class of Bethe 
Ansatz equations that contained the Richardson solution and the Lipkin-Meshkov-Glick 
(LMG) model. 

In recent years, many papers have appeared devoted to finding new examples of BCS 
systems solvable by the Bethe Ansatz [4T14"5"9"] . In contrast to the Richardson solution for 
s-wave pairing, several of these newer systems exhibit quantum phase transitions which 
can be identified by a change in the character of the Bethe roots as a coupling constant 
is varied. Such behaviour has also recently been observed in a different context of 
bosonic models [60J. The purpose of the present paper is to provide a first step towards 
extending the application of numerical methods based on the BA/ODE correspondence 
to these more recent BCS systems. In particular we will conduct an analysis for the 
degenerate two- level limit. 

The BCS Hamiltonians we consider below are generally expressed in terms of 
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fermions. Given the blocking effect |29J of the pairing interaction, for simplicity we 
hereafter only consider for all models the subspace of unblocked states. This effectively 
reduces the dimension of the Hilbert space of states for the model. In this case it is more 
convenient to write the Hamiltonians only in term of hard-core boson operators. The 
hard-core boson operators &j and bk (Cooper pair creation and annihilation operators) 
satisfy the following commutation relations 



we further effectively reduce the dimensionality of the Hilbert space by only considering 
certain "unblocked symmetric states" which will be defined later. One motivation for us 
to do this is that we want to illustrate that the numerical methods employed here provide 
solutions of the Bethe Ansatz equations for excited states as well as the ground state. 
In order to display this data in a transparent manner, it is an advantage to restrict our 
analysis to a subspace of the full Hilbert space. If the ground state on the unrestricted 
Hilbert space is non-degenerate, symmetry arguments lead to the conclusion that the 
ground state is an unblocked symmetric state. Finally, we mention that the method 
employed here (and in (43j,[44] ) is distinct from other approaches [32 |l36ll39t l42 | [53 t l54] 
in that we do not track solutions of the Bethe Ansatz equations from the zero coupling 
limit. In this respect issues surrounding critical points are avoided. 

In Section 2 we will introduce a class of nonlinear algebraic equations and obtain the 
corresponding ordinary differential equations involving only polynomials. We classify the 
different cases and discuss the application to two-level models. In Section 3, we discuss 
an algorithm to numerically generate the corresponding generalised Heine-Stieltjes and 
Van Vleck polynomials. We apply this approach to four systems: the Richardson s- 
wave model [25], the p + ip-wave pairing model [SllESllHI], the p + ip-w&ve pairing 
model coupled to a bosonic molecular pair degree of freedom [56], and a new solution 
for a d + id-wave pairing model with additional interactions [EI] . We will compare the 
ground-state roots obtained by this numerical method with the theoretical curve for an 
arc obtained in the continuum limit by a singular integral equation. In two cases we 
observe a phenomenum whereby all the ground-state roots collapse at the origin at a 
particular value of the coupling parameter. The method also allows us to obtain the 
roots for the excited states, which is important for studies of the dynamics analogous 
to those in [3T|[38] . 

2. The generalised Heine-Stieltjes correspondence 

The systems studied previously using the polynomial approach [3l ^H2TPH] belong to the 
class of coupled nonlinear algebraic (Bethe Ansatz) equations of the form 






(2) 
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where C,D,pi and £j are real parameters. The parameters can be interpreted in the 
context of BCS systems as single particle energy levels, L is the number of levels and 
M is the total number of Cooper pairs. The various systems studied in [50 |l53l - l56| [59] 
have Bethe Ansatz equations that belong to a different class than the one given by Eq. 
02]), taking the following form 

L M . 

V^^-2V^— + 4 + - + C = 0, Z = 1,...,M. (3) 

j^yi-et j^yi-yj yf vi 

In the equation above, the real parameters A, B and C may depend on L and M or 
other constants depending of the nature of the physical problem considered. 
We start by constructing the following polynomials 

M L 

Q(z) = H(z- yj ), P(z) = H(z-e i ). (4) 

j=l i=l 

It is well known that such polynomials satisfy the following relations 

Q'ivi) j^vi-v,' p (yi) t(yi- £ i' 

We introduce a polynomial W{z) which is of order L — 1 such that for the set of real 
parameters pi and £j it satisfies 

W(z) J~, p, 



E 



P(Z) Z-Ei 

and in the case pi = 1 Vi, as can be seen from Eq.fjSJ), W{z) = P'{z). These relations 
can be used to construct from Eq. (jHJ) the following differential equation 

A B W( yi ) Q"( yi ) _ 

yf m P{yi) Q'{yi) 

This equation can be written in the form A 2 (yi)Q"(yi) + Ai(yi)Q'(yi) = 0. Because, the 
polynomial Q(z) vanishes at the solutions y\ of the Bethe Ansatz equations we can thus 
form the following second order differential equation 

A 2 {z)Q"{z) + A l {z)Q'{z) = A Q (z)Q(z), (6) 

where the order \Aj\ of the polynomial Aj(z) depends on when the parameters A, B, C 
vanish. Seven cases can occur, as given in Table 1. 

For a given L, the polynomial A (z) can be seen as a generalised Van Vleck 
polynomial and the polynomial Q(z) as a generalised Heine-Stieltjes polynomial. For our 
numerical investigations we will study degenerate two-level models, where each of the 
distinct levels E\ and e 2 have have degeneracy L/2. However the method can be applied 
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Table 1: Cases of generalised Heine-Stieltjes and Van Vleck polynomials for the 
differential equation given by Eq. dSJ) 



Case 


A 


B 


C 


A 2 


Ax 




\A 2 




\Ai 




\Ao\ 




1 


A^O 


B^O 


C^O 


z 2 P 


(-A- 


Bz-Cz 2 )P-z 2 W 


L + 


2 


L + 


2 


L + 


1 


2 


A^Q 


B^O 


C = 


z 2 P 


(-A- 


Bz)P-z 2 W 


L + 


2 


L + 


1 


L 




3 


A^Q 


B = 


C ^ 


z 2 P 


(-A- 


Cz 2 )P-z 2 W 


L + 


2 


L + 


2 


L + 


1 


4 


A^Q 


B = 


c = o 


z 2 P 


-AP- 


- z 2 W 


L + 


2 


L + 


1 


L 




5 


A = 


B^O 




zP 


(-B- 


Cz)P-zW 


L + 


1 


L + 


1 


L 




G 


A = 


B^O 


C* = 


zP 


-BP- 


- zW 


L + 


1 


L 




L - 


1 


7 


A = 


B = 




P 


-CP- 


- W 


L 




L 




L - 


1 



to general L-level models at the cost of much more involved numerical calculations. For 
our study the problem becomes mathematically equivalent to taking p« = L/2 Vi for the 
first term of Eq. fl3]). The Bethe Ansatz equations then take the following form 

A B „ L ( 1 1 \ A 1 

yf vi 2\£i-vi s 2 -yij j^yi-yj 

3. Numerical method based on B A/ODE correspondence and examples 

We will use Eq. ([6]) to construct the polynomials Q(z) and A (z) by inserting expansions 
of these polynomials and solving the corresponding system of equations. From the 
coefficients of the polynomial Q(z) we obtain the roots by standard techniques, and 
they correspond to the solutions of the corresponding Bethe Ansatz equations. This 
approach is similar to the one taken in 

We start by taking the following expansions (where M and L are fixed) 

M K 
3=0 3=0 

with ajvf = 1 and the remaining otj and /3j are coefficients to be determined numerically. 
Inserting the expansions given by Eq. (jHJ) into Eq. ([6]) yields two matrix equations that 
we need to solve. The coefficients of z % for i = 0, M generate a (M + 1) x (M + 1) 
matrix F that satisfies the equation Fv = /3 v with v r = (a , •••,«m)- The coefficients 
of z % for % — M + 1, ...,M + K will generate a K x (M + 1) upper triangular matrix 
P that satisfies Pv = 0. The matrix entries are all linear in the coefficients of the 
generalised Van Vleck polynomials {3k}- We obtain a set of coefficients otk (with 

k G [0, M]) for each solution set of the Bethe Ansatz equations and a corresponding set 
of coefficients (3j (with j e Some of the (3j may be independent of the and 

thus do not depend on the set of Bethe roots. We have adapted a Mathematica code 
discussed in [43J to implement the calculations. 
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An effect of restricting to the unbloacked symmetric states of the degenerate, two- 
level models which are governed by Eq. (j7]) is that the order of the polynomials A Q , A\ 
and A2 is fixed as we change the value of L and M. The polynomials are respectively 
of order 3, 4 and 4 at most, and thus the parameter K is at most 3: 

4 4 

j=0 j=0 

When A\ reduce to a polynomial of order 3, the equation fall into the class of differential 
equations recently studied in [17]. Once the c*i and (3j are obtained, each solution 
corresponding to solution set of Bethe roots, the next step is to compute the roots from 
the generalised Heine-Stieltjes polynomials. This can be done in principle using one of 
the many well-known methods implemented in available softwares, and for an arbitrary 
order. But in practice this is a more complicated issue. A difficulty that appears is 
that these generalised Heine-Stieltjes polynomials will be polynomials of order M and 
care is needed as it is known that numerical methods to find roots of polynomials can 
have instabilities [621464] . This problem can even affect the structure of the roots (i.e. 
real vs. complex conjugate pairs). We will adopt in this paper an approach based on 
monomial expansion but we will take care of using a sufficient working precision. We 
will show that the method appears to be reliable and robust. We will present explicit 
examples of polynomials in the Appendix A, illustrating that the coefficients can be of 
very different orders of magnitude. 



3.1. Example 1: Richardson s-wave pairing model 

We first consider the exactly solvable s-wave pairing BCS model [25H36] given by the 
following Hamiltonian expressed in terms of the hard-core boson operators 

L L 

H = ^jN J -GY J b]b k . 

J'=l j,k 

Omitting blocked states, this Hamiltonian acts on a Hilbert space of dimension 2 L . The 
energy is given by 

M 

where the roots yj satisfy Eq. fl3]) with C = G^ 1 and A — B — 0. This corresponds to 
case 7 of Table 1. 

Imposing L is even the degenerate two-level limit is obtained by setting 



e x j odd, 
e 2 j even 
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in which case the Hamiltonian may be expressed as 

H = L + 2e x Sl + 2e 2 S 2 2 - G{S+ + <S+)(Sr + s 2 ) 

where 



si 


j odd 


s 2 


3 even 


(10) 


si 


= E«*. 

j odd 


s 2 


= E«*. 

3 even 


(11) 


st 


= E 6 J. 

3 odd 


st 


= E»J 

3 even 


(12) 



provide two representations of the su(2) algebra. Now the energies given by (Q, where 
the roots yj satisfy Eq. ((Tj) with C = G~ x and D = 0, only provides a subset of 
the spectrum for the Hilbert space of dimension 2 L . These correspond to states which 
are invariant under the mutual interchange of the even subscripts, and the mutual 
interchange of the odd subscripts, which label the components of the tensor products. 
We term these unblocked symmetric states. This subspace has dimension (L/2 + l) 2 , 
and will be the focus of our study in this section. In representation-theoretic terms, the 
hard-core operators satisfying Eq. (CQ) can be mapped to L irreducible spin- 1/2 su{2) 
representations, whereas (I10|I11|I12I) are two reducible su(2) representations. The Eq. 
(J7j) with C = CE 1 and D = is only associated with one irreducible component of these 
reducible representations, that with spin-L/4 for each su{2) copy. Similar considerations 
will apply for all subsequent models that we study. 

Setting 7 = ei + e 2 and r] = S\S 2 , in the degenerate two- level limit the differential 
equation takes the form 

(z 2 - 1Z + V )Q" + (-±z 2 + (l- L )z+^yL- Q' - {Piz + fo)Q = 0. 

From the first and second terms in the expansion of the polynomial Q (using Eqs. (H|9|) ) 

M 

Q(z) = z M - Ez M ~ l + ... + (-1) A/ y j} 

i=i 

we have that 

/3i = ~, P = - (^M +^ + M + LM-M 2 y 

Because the parameter (3 depends on the energy, we obtain a Van Vleck polynomial 
corresponding to the Heine-Stieltjes polynomial of each eigenstate of the Hamiltonian. 

We take the case where £2 = 1 and S\ = —1, following j27J[32]. We choose L = 100, 
a system at half-filling M = 50, and introduce the scaled coupling constant g = GL 
which will be used throughout. From the numerical method based on the BA/ODE 
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Re(;y) Re(y) Re(y) 



(a) (b) (c) 

Figure 1: Roots for the ground state of the degenerate, two-level Richardson model with 
e 2 = \ )£l = -1, L = 100, and M = 50: (a) g = 1/2, (b) g = 1, and (c) g = 3/2. Also 
shown are the theoretical curves derived in Appendix B for the large- L limit. 




Figure 2: Roots for all unblocked symmetric states of the degenerate, two- level s-wave 
pairing model with e 2 = 1, e± = —1, L = 100, and M = 50: (a) g = 1/2, (b) g = 1, 
and (c) g = 3/2. In each case a total of 51 sets of roots are displayed, where each set 
contains 50 roots. 



correspondence we reproduce the ground state results of [271132] . We see that the results 
obtained by this method also agree with the theoretical distribution curve in the large- L 
limit which is calculated in Appendix B. The case g = 1 is identified as a critical point 
at which the character of the arc changes from being a closed curve to an open curve. 
The method we use also computes the roots of the excited states. Fig. 2 shows all sets 
of the Bethe roots for given values of g. The number of states in the M = 50 sector is 
51, so the number of data points in each panel is 50 x 51 = 2550. 
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Table 2: Phase diagram of the p + ip-w&ve pairing model 





Phase 


Constraint between coupling constant g and filling fraction x 


1 


Weak coupling BCS 


x > 1 - g- 1 


2 


Moore-Read line 


xmr = 1 - 


3 


Weak pairing 


(1 - g~ x )/2 < x < 1 - g- 1 


4 


Read- Green line 


xrg = (l - g- 1 )/^ 


5 


Strong pairing 


x<(l-g- 1 )/2 



3.2. Example 2: The p + ip-wave pairing Hamiltonian 

Next we examine the p + ip-wave pairing model, whose exact solution has only recently 
been derived [50~|l53]. The connection to Heine-Stieltjes and Van Vleck polynomials for 
this model was noted in [51]. Up to a canonical transformation, the Hamiltonian reads 



L L 

i=l j<k 

A notable feaure of this model is the ground-state phase diagram, which is summarised 
in Fig. 2. The phase boundaries known as the Moore- Read line and the Read-Green 
line have the property that they are independent of the parameters Sj. Consequently, 
signatures of the phase boundaries should still be present in the degenerate two-level 
limit. 

The energy of each eigenstate is given as the sum 

M 

E=(l + 0)^2 yj , (13) 

where the roots yj satisfy Bethe Ansatz equations of the form given by Eq. ([3]) with 
B = G- 1 - L + 2M - 1 and A = C = 0. It thus belongs to case 6 of Table 1. In the 
degenerate two-level limit, the corresponding second order differential equation is 

(Z 3 _ 1Z > + nz)V + A,(-l + G(l + ^-2M)) + 7(2 - G (2 + L-4M)), 

y G 2G 



G 

Using the expansion of the polynomial Q(z) and the energy expression (flBl leads to 

M , a ( E 1 M , 1 LM ,,2 



We consider the case e 2 = 1, £\ = 1/2, L — 200, and M = 50. Ground state roots 
are depicted in Fig. 3, along with the theoretical curves derived in Appendix B for the 
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Figure 3: Roots for the ground state of the degenerate, two-level p + ip-w&ve pairing 
model with e 2 = 1 and E\ — 1/2, L — 200, and M = 50: (a) g — 1/2 (Weak coupling 
BCS), (b) g = 4/3 (Moore-Read line), (c) g = 3/2 (Weak pairing), and (d) g = 2 
(Read-Green line). Also shown are the theoretical curves derived in Appendix B for the 
large- L limit. 



large-L limit. A curious feature here is case (b) corresponding to the Moore-Read line. 
At this point there is a change from an open curve to a closed curve in the continuum 
limit, however all the roots collapse to the origin in any finite system. This is a result 
of the limits g — > qmr and L — > oo not commuting, and was identified in [5QII53] as a 
zeroth order quantum phase transition. This is in stark contrast to Fig. 1 (b) . We also 
highlight that for case (c) the curve consists of two connected components, one of which 
is a closed curve and another which is open. At the Read-Green line shown in (d) the 
closed curve has contracted to a point at the origin. The theoretical curves shown in 
Fig. 3 have exactly the same qualitative features as is found in the general p + zp-wave 
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Figure 4: Roots for all unblocked symmetric states of the degenerate, two-level p + ip- 
wave pairing model with e 2 = 1 and e\ = 1/2, L = 200, and M = 50: (a) g = 1/2, 
(b) g = 4/3, (c) g = 3/2, and (d) g = 2. In each case a total of 51 sets of roots are 
displayed, where each set contains 50 roots. 



pairing model [53, 53]- ^ s ^ n the previous example, the number of states in the M = 50 
sector is 51, so the number of data points in each panel of Fig. 4 is 2550. 

3.3. Example 3: The p + ip-wave pairing Hamiltonian coupled to a bosonic 
molecular pair degree of freedom 

We next consider the case of the p + ip-wave pairing Hamiltonian coupled to bosonic 
molecular pair [56]. This model is given by the following Hamiltonian 

L L L 

H = J2 e i N i - F 2 GNq V^{b)b k + blbj) \fc(bob] + b\b 3 ) 

j=l j<k j=l 
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where bo, b\ satisfy the bosonic commutation relation 

K b Q ] = I 

and No = b^bo- This system has the energy spectrum 

M 

E = (l + 0)^2 yj , 

where the roots yj satisfy the Eq.© with C = 0, B = G' 1 + 2M — L — l and A — F 2 , 
thus belonging to case 2. We form the following differential equation 

(z 4 - -fz 3 + r]z 2 )Q" + (-F 2 ri + (F 2 7 + r] - £ + r]L - 2i]M)z 

G 

+ (-F 2 - 7 + 1 - ^ + 2 1 M)z 2 + (1 - i - 2M)z 3 )Q' - (f3 2 z 2 + (3 lZ + (3 )Q = 0. 

The coefficient of the generalised Van Vleck polynomial can be rewritten in the following 
form 

h = -„a +U \, ft- f2E + (-2( 1 -F*G) +1 GL)M-2 7 GM-- 



G J V 2G 

In this instance, some coefficients of A\ and all coefficients of A depend of M. The 
coefficient /3o is more complicated than for the previous examples, as it involves a double 
sum of the roots and is not given only in terms of the energy. 

We consider the choice e 2 = 1 and e t = 1/2, L = 32, M = L/2 and F = Vl28 to 
undertake the numerical calculations. Ground-state roots are shown in Fig. 5. They 
always lie on the negative real-axis. The roots of all unblocked symmetric states are 
shown in Fig. 6. There are no significant qualitative changes in the pattern of roots as 
the coupling parameter g is varied. 

3.4- Example 4 : An extended d + id-pairing Hamiltonian 

Using the results concerning the conserved operators of the Richardson model, it is 
possible through a change of variables for the single particle levels to obtain an integrable 
Hamiltonian given by [61] 



L L 

H = CM - G eMfyk + h A + 1N 3 N k ). 

i=l j,k=l 

In the absence of the NjN^. interaction terms, this model maps through a canonical 
transformation to a d + id-pairing Hamiltonian. 
The energies are given by 

M MM L 

E = Y,m-2Gj2J2y^- G Y. £ r 

1=1 1=1 tfl j=l 



Heine- Stieltjes and Van Vleck polynomials associated with integrable BCS models 13 



-0.6 -0.4 
Re(y) 



-0.2 _ 

-20 



(a) 



-10 

Re(y) 



(b) 



-1000 -800 -600 -400 -200 

Re(v) 



(c) 



Figure 5: Roots for the ground state of the degenerate, two-level p + ip-w&ve pairing 
model coupled to a bosonic molecular pair with e 2 = 1, £\ — 1/2, L = 32, M = 16, and 
F = V128: (a) g = 1/10, (b) g = 1, and (c) g = 10. Also shown are the theoretical 
curves derived in Appendix B for the large- L limit. 



0.03 
0.02 
0.01 
; 0.00 
-0.01 
-0.02 
-0.03 



-0.5 



0.0 0.5 

Re(y) 



(a) 



0.05 



().()() 



-0.05 



o.ior 



0.05 



0.00- 



-0.05 



-0.10L 



1.0 



-2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 
Re(y) 

(b) 



-2 -1 

Re(j) 



(c) 



Figure 6: Roots for all unblocked symmetric states of the degenerate, two-level p + ip- 
waved pairing model coupled to a bosonic molecular pair with e 2 — 1, E\ — 1/2, L = 32, 
M = 16, and F = Vl28: (a) g = 1/10, (b) g = 1, (c) g = 10. In each case a total of 
153 sets of roots are displayed, where each set contains 16 roots. 



where the set of parameters yi satisfy Eq.()3]) with C — 0, B — 2M — 2 — L and 



L 



case thus belong to the previous example. We will 



see however, that the parameters A, B and C will generate very different behaviours for 
the roots of the generalised Heine-Stieltjes polynomials. The equation can be put into 
the form 

(, 4 - 7 z 3 + V z*)Q" + (- JL + l - ir} L + Un + ^ - (7 ' + 2 2V)L - 2 VM ) z 



+ -27 



^ + 2 1 M^j z 2 + (2- 2M)z^ Q' - (f3 2 z 2 + ftz + f3 )Q = 0. 



Heine- 
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— x — x— X— x-> 



x x x x x 
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Re(y) 
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-0.5 0.0 
Re(j) 



(b) 



0.03 " 
0.02 
0.01 
0.00 
-0.01 
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-0.03 _ 
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Re(y) 



(c) 



Figure 7: Roots for the ground state of the degenerate, two-level extended d + id-wave 
pairing model with e 2 = 1, e x = 1/2, L = 64, and M = 32: (a) g = 49/75 , (b) g = 2/3, 
and (c) g = 51/75. Also shown are the theoretical curves derived in Appendix B for the 
large- L limit. In case (b) the curve has contracted to a point at the origin. 



We can show that the coefficients of the Van Vleck polynomial may be written as 
/3 2 = -((M-l)M), /3 1 = -( 7 M+^- 7 M 2 ), 

Po = ~ ( 2! + ¥ " VM ~ 2^" + \ (7 ~ v) LM + VM2 ) ■ 

For the numerical calculations we choose e 2 — 1, £\ — 1/2, L = 64, M = 32. The 
total number of states is given by M. Results are shown in Figs. 7-9. These figures 
illustrate the changing behaviour of the roots near the point g = 2/3. The the ground- 
state roots form an arc in the regions near this point and collapse at the origin at the 
this point (similar behavior as that of the Moore- Read line of the Example 3). 

We finally identify four ground-state phase boundary lines which are associated 
with changes in the topology of the root distribution curve in the continuum limit. 
These cases are respectively shown by Fig. 1(b) for the s-wave model, Fig. 3(b) and 
Fig. 3(d) for the p + ip-wave model, and Fig. 7(b) for the extended d + id-wave model. 
From the the numerical calculation using the BA/ODE correspondence, we find that 
the behaviour of these curves and the ground-state roots of the ground state can differ 
greatly. Remarkably, each of the cases show distinctive features which suggests that the 
four boundary lines should exhibit intrinsically different consequences for the properties 
of the models. The observations are summarised in Table 3. Another curious observation 
concerns the excited states. At the Read-Green line for the p+ip-paihng Hamiltonian as 
depicted in Fig. 4(d), and at the boundary for the extended d + id-pairing Hamiltonian 
as depicted in Fig. 9(b), the roots for all excited states cluster near common curves. 
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0.01 



-0.01 




0.00 

Re(.y) 



(a) 



0.02 



0.01 



0.00 



-0.01 



-0.02 



-0.02 




-0.01 



0.00 



0.01 



0.02 



(b) 



Figure 8: Theoretical curves derived in Appendix B for the large-L limit of the 
degenerate, two-level extended d + id-wave pairing model with £ 2 = 1, £i = 1/2, and 
x = lim M/L = 1/2: (a) From left to right g = 203/300, g = 202/300, # = 201/300, 

L— >oo 

(b) From left to right g = 199/300, g = 198/300, g = 197/300. The limiting behaviour 
indicates that the curve contracts to a point at the origin when g = 2/3. 
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Figure 9: Roots for all unblocked symmetric states of the degenerate, two-level extended 
d + id-wave pairing model with e 2 — 1, £1 = 1/2, L = 64, and M = 32: (a) g = 49/75, 
(b) g = 2/3, and (c) g = 51/75. In each case a total of 33 sets of roots are displayed, 
where each set contains 32 roots. 

Table 3: Behaviour of the ground-state roots in a finite-sized system, and the distribution 
curves in the continuum limit, for ground-state phase boundary lines 



Model 


Figure 


# roots at the origin 


Behavior of the curve 


s 


1(b) 


Nil 


Closed/open transition 


p + ip (Moore- Read line) 


3(b) 


All 


Open/closed transition 


p + ip (Read-Green line) 


3(d) 


Nil 


Closed curve collapse 


Extended d + id 


7(b) 


All 


Open curve collapse/revival 
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4. Conclusion 

We used the BA/ODE correspondence and the related generalised Heine-Stieltjes and 
Van Vleck polynomials to numerically obtain the roots of the Bethe Ansatz in the 
limit of two levels for four cases of integrable pairing models. We reproduced results 
obtain by other methods concerning the ground states, and we were able to study the 
behavior of the roots for the excited states as we change the coupling constant. The 
CPU time needed was consistent with what was mentioned in [43J. In this paper we 
used a monomial expansion and found that the method is robust and reliable in all 
cases with a moderate value of M and a sufficient working precision. The application 
to general A-level models is straightforward, however, the CPU time needed to perform 
the calculations will increase. The most efficient method to obtain the roots of these 
polynomials in such cases needs further study. The existence of recursion operators as 
suggested in [H] is also an aspect to be investigated as it could lead to faster methods. 
This approach is not limited to integrable BCS models but can be applied to integrable 
systems that appear in other contexts, and to other classes of BA equations. 

An interesting aspect to investigate is the relation between the Schr'odinger form 
of the equations (jSJ) for these four cases, and quasi-exactly solvable systems (QES) [65J. 
The classification of QES systems was performed for certain classes of equations of 
the form given by (J6j) with polynomials of order 2,3 and 4. The equations in this 
paper seem to be beyond this classification, and also other systems obtained more 
recently [65T - I67] . Supersymmetric quantum mechanics appears to be useful method 
to obtain such systems [6"6|l6"T]. 
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Appendix A. Explicit generalised Heine-Steiltjes and Van Vleck 
polynomials 

The following cases correspond to the four Examples of Sect. 3. It is seen that the 
coefficients can differ across many orders of magnitude. The numerical precision used 
in the computation was significantly higher than the number of digits shown below. 
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Example 1: ground state at g — 1 



a 


= 4.3 






Oti 


— 


2.1 


X 


10 2 , 


a 2 




4.9 


X 


10 3 , 


CK3 


— 


7.3 


X 


10 4 , 


a>4 


= 8 


1 


X 


10 6 , 


a 5 


— 


7.1 


X 


10 6 , 




— 


5.1 


X 


10 7 , 




— 


3.1 


X 


10 8 , 


a$ 


= 1 


6 


X 


10 9 , 


a 9 


— 


7.1 


X 


10 9 , 


«10 


— 


2.8 


X 


10 10 , 


a n 


— 


9.8 


X 


10 10 , 


Ot-Yl 


= 3 


1 


X 


10 11 , 


"13 


— 


8.6 


X 


10 11 , 


"14 


— 


2.2 


X 


10 12 , 


"15 


— 


5.1 


X 


10 12 , 


«16 


= 1 


1 


X 


10 13 , 


Ctyj 


— 


2.1 


X 


10 13 , 


«18 


— 


3.7 


X 


10 13 , 


«19 


— 


6.0 


X 


10 13 , 


«20 


= 9 


1 


X 


10 13 , 


OLl\ 


— 


1.3 


X 


10 14 , 


«22 


— 


1.6 


X 


10 14 , 


«23 


— 


1.9 


X 


10 14 , 


«24 


= 2 


1 


X 


10 14 , 


"25 




2.1 


X 


10 14 , 


«26 


— 


2.0 


X 


10 14 , 


«27 




1.7 


X 


10 14 , 


«28 


= 1 


4 


X 


10 14 , 


«29 




1.0 


X 


10 14 , 


"30 




7.0 


X 


10 13 , 


«31 




4.4 


X 


10 13 , 


«32 


= 2 


6 


X 


10 13 , 


«33 




1.4 


X 


10 13 , 


"34 




6.7 


X 


10 12 , 


"35 




3.0 


X 


10 12 , 


"36 


= 1 


2 


X 


10 12 , 


a 37 




4.5 


X 


10 11 , 


"38 




1.5 


X 


10 11 , 


"39 




4.6 


X 


10 10 , 


"40 


= 1 


2 


X 


10 10 , 


«41 




2.9 


X 


10 9 , 


«42 




6.2 


X 


10 8 , 


«43 




1.1 


X 


10 8 , 


0:44 


= 1 


8 


X 


10 7 , 






2.3 


X 


10 6 , 


«46 




2.5 


X 


10 5 , 






2.1 


X 


10 4 , 












«48 




1.3 


X 


10 3 , 


«49 




5.1 


X 


10 4 , 













(3o = 2.5 x 10 3 , ft = -5000. 



The coefficient (5\ is an exact value which is independent of the state. 



Example 2: ground state at g = 3/2. 



a = 


0, 










= 0, 




a 2 




0, 




«3 


= 0, 








«4 = 


0, 










= 0, 








0, 




0:7 


= 0, 








a 8 = 


0, 










= 0, 




«10 




0, 




a n 


= 0, 








«12 = 


-1 


7 


x 10" 


-46 


«13 


= -6.5 x 10" 45 , 


«14 




-2.2 x 10" 


43 


"15 


= 7.3 


X 


10" 


-42 


a 16 = 


-2 


3 


x 10" 


-40 


an 


= 6.9 x 


IO- 39 , 


«18 




-2.0 x 10" 


37 


«19 


= 5.8 


X 


10" 


-36 


«20 = 


-1 


6 


x 10" 


-34 


«21 


= 4.6 x 


10" 33 , 


«22 




-1.3 x 10" 


31 


«23 


= 3.5 


X 


10" 


-30 


«24 = 


-1 





x 10" 


-28 


"25 


= 2.9 x 


10- 27 , 


«26 




-8.9 x 10" 


26 


«27 


= 2.9 


X 


10" 


-24 


«28 = 


-1 





x 10" 


-22 


«29 


= 4.1 x 


10- 21 , 


"30 




-1.9 x 10" 


19 


«31 


= 1.2 


X 


10" 


-17 


«32 = 


-1 


2 


x 10" 


-15 


"33 


= 4.3 x 


10- 13 , 


"34 




7.1 x 10 -11 , 


"35 


= 4.3 


X 


10" 


-9 

i 


"36 = 


1.4 


X 


10- 7 




"37 


= 3.2 x 


10" 6 , 


"38 




5.0 x 10~ 5 , 




"39 


= 5.7 


X 


10" 


-4 

j 


"40 = 


5.0 


X 


10~ 3 




a 41 


= 3.3 x 


10- 2 , 


«42 




1.7 x 10 _1 , 




"43 


= 6.6 


X 


10" 


-1 

) 


0:44 = 


2.0 








«45 
«48 


= 4.7, 
= 8.3, 




«46 
«49 




8.1, 
4.3, 




«47 


= 1.0 


X 


10 1 





Po = 6.8 x 10 4 , 



/3i = -27500/3. 
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The coefficient 0i is an exact value which is independent of the state. 



Example 3: ground state at g = 1. 



a = 2.4 x 10 
«4 = 0.0, 
a 8 = 0.23, 
"l2 = 11, 



-ii 



Oil 
"13 



2.3 x HT y , 
6.0 x 10~ 4 , 

9.4 x 10 -1 , 
13, 



1-7 

o 6 = 6.1 x 10~ 3 , 

"10 



a 2 = 1.1 x 10" 



2.9, 



Ol4 = 10, 



"3 
O'v 

on 

"15 



2.8 x 10" 
0.043, 
6.6, 
4.8, 



Po = -2.0 x 10 2 , 0t = 8.6 x 10 2 , fa = -768. 

The coefficient 02 is an exact value which is independent of the state. 



:ampl 


e 4: ground state at g = 


= 102/150. 


















"o = 


6.1 x 10~ 55 , 


ax — 6-5 x 


io- 52 , 


"2 


= 3.5 x 


io- 


-49 


"3 


= 1.3 


X 


10" 


-46 


«4 = 


3.3 x 10" 44 , 


« 5 = 7.1 X 


lO" 42 , 


"6 


= 1.3 x 


10" 


-39 


"7 


= 1.9 


X 


10" 


-37 


"8 = 


2.4 x 10~ 35 , 


a 9 = 2.7 x 


10" 33 , 


"10 


= 2.7 x 


10" 


-31 


"11 


= 2.4 


X 


10" 


-29 


"12 = 


1.9 x 10~ 27 , 


«13 = 1.4 X 


IO" 25 , 


"14 


= 9.0 x 


10" 


-24 


"15 


= 5.3 


X 


10" 


-22 


"16 = 


2.8 x 10~ 20 , 


Oi7 = 1.4 X 


io- 18 , 


"18 


= 5.9 x 


10" 


-17 


"19 


= 2.3 


X 


10" 


-15 


"20 = 


8.2 x 10~ 14 , 


o 2 i = 2.6 x 


io- 12 , 


"22 


= 7.3 x 


10" 


-11 


"23 


= 1.8 


X 


10" 


-9 


"24 = 


4.0 x 10~ 8 , 


"25 = 7.7 X 


io- 7 , 


"26 


= 0.0, 






"27 


= 1.0 


X 


10" 


-4 

; 


"28 = 


1.9 x 10~ 3 , 


a 29 = 1.7 x 


10" 2 , 


"30 


= 1.1 X 


10" 


-1 

5 


"31 


= 4.6 


X 


10" 


-1 

) 



O = 5.1 x 10 2 , X = -304/17, 02 = -1-0 x 10 3 . 

The coefficients 0i and 02 are both exact values which are independent of the state. 

Appendix B. Solution of the Bethe Ansatz equations in the continuum limit 

If in the limit L — > oo the roots of the Bethe Ansatz equations ([3]) are densely distributed 
on a curve T in the complex plane, we can look to find a solution via integral equation 
methods. We refer to this as the continuum limit. Specifically we obtain the singular 
integral equation 



A B r p{e) 

y 2 y Jn y-e 



,Mv') 



o 



y 



(B.l) 



where Q denotes the interval of the real line where the energy levels e lie, distributed 
according to a density p(s) such that 



de p(e) 



1. 
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These integral equations may be solved using complex analysis techniques (cf. 
|33],|53j|5ll|6T]). Since the energy must be real, we assume that the arc Y is invariant 
under reflections about the real axis. The solution for V having end points a = e — i5 
and b = e + i5 is of the form 

r(y)\dy\ = ^-.(h + (y) - h_(y))dy, 



de 



h u H h -r 

e-y y y 2 



Ky) = R(y) 
R(y) = y/(y- a )(y-b) = ^{y-ef + 5 2 . 
The Cauchy principal value in Eq. ( IB. II) can be written as: 

„ r(y') f dy' h(y') 



Jt V - V Jc r y'-y 

where Cr is a closed curve which contains T in the interior. Similar equations apply for 
the case where T lies on the real axis with support on the interval [a, b}. 

The solution can thus be obtained on the form of contraints that involve integrals 
of the density p{e). For all four cases relations of the form 



T x (a, b), 5f = J r 2 (a,&) 



(B.2) 



can be obtained. For the two-level models we use p(e) = 5(e — ei)/2 + 5(e — e 2 )/2, and 
for a given g and x we may numerically determine the end points of T. From these we 
use the equation 



Re 







to obtain the arc. 

Alternatively, the arc may form a closed curved in the complex plane. In this case 
the solution of the singular integral equation takes the form 



s(y) 

and the curve is given by 



1 

m 



ds 



P (g) , , V l , V 2 

V v H V — 

e-y y y 2 



Im 



o. 



Here, the introduction of a integration constant is necessary. It can be obtained by 
imposing the constraint that the value w such s(w) = lies on the curve. In the case 
that T = Ti U T 2 , where Ti is a closed curve that touches a point b of T 2 = (a, b) ^ Q 
we have 



s(v) 



V(y ~ a )(y- h ) 



in 



de h v H h 

e-y y y 2 
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Below we provide details for when the arc V is associated to the ground-state roots 
in each of the four examples. For each case it is found that (fi(e) = p(e)/\R(e)\. 

Example 1 

The corresponding equation for Eq. (1B.2I) are 

« = ^<M«)(i-^), 



9 Jn R(e) 

with uo = u\ = ui = 0. In the case of closed curve we have v 2 = v 1 = and vq = l/(2g). 
In this case explicit expressions for the curve can be obtained for the open arc 



2 1 2 1 2 ^ x 9 v. 

X +V +9 = tanh(2^) ' 



and for the closed arc 



990 2^t/<s^ 
X +V +£ ' = ta, lh (2(x- & ) 9 -.)' 9<£l 



with 



c 9 , / £1 + Xo \ L 9 

2 \£i-XoJ V £ i 



Example 2 

The corresonding equations for ( IB. 21) are 



2x = 1 - - + Ve 2 + S 2 [ de 
9 Jn 

ep(e) 



R(s) 



de ■ 



9 Jn R{ £ ) 

with m = «2 = and u\ = (2s — 1 + g^ 1 )(e 2 + 5 2 ) 3 ^ 2 /2. In the case of a closed 
curve we have v = v% — and t>i = 2a; — 1 + g -1 . In the case of a closed curve I\ 
that touches a point 6 of T 2 = (a, 6) ^ we have v = v 2 = and t>i = (2a; — l + g _1 )/a&. 

Example 3 

In this case the arc associated with the ground-state roots is always an interval [a, b] on 
the real axis. The corresponding equations for (1B.2|) are 



2^G~b Jn 2R(e) 

1 P £ P( £ ) 



g yfal, J n 2R(e) 



+ / de 
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with 

q = 2x-l + -, 
9 

F 



f y/V 

Uq = 0, 

1 / f 2 (a + 6) 
Mi = = q + 



u 2 



P 



2Vab 



Example 4 

In this case the equations given by ( IB. 21) take the form 



I puj red 

- = 2 deep(e) + 2Ve 2 + 5 2 / deeMe) 
9 Jo Jo 



with 



c = 2x — 1, 

f =—- I deep{e), 
2g Jn 

U = 0, 

Ml = + 



'U 2 



V?TF (e 2 + 5 2 ) 3 / 2: 
/ 

Ve 2 + 5 2 ' 
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